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Abstract 

The popular QR algorithm for solving all eigenvalues of an unsymmetric matrix is 
reviewed. Among the basic components in the QR algorithm, it has been concluded from 
this study, that the reduction of an unsymmetric matrix to a Hessenberg form (before 
applying the QR algorithm itself) can be done effectively by exploiting the vector speed and 
multiple processors offered by modern high-performance computers. 

Numerical examples of several test cases have indicated that the proposed parallel-vector 
algorithm for converting a given unsymmetric matrix to a Hessenberg form offers 
computational advantages over the existing algorithm. The time saving obtained by the 
proposed method is increased as the problem size increased. 


I. Introduction 


The algorithms for symmetric matrices [1-3] are highly satisfactory in practice. By 
contrast, it is impossible to design equally satisfactory algorithms for the nonsymmetric cases, 
which is needed in Controls-Structures Interaction (CSI) applications [1,4]. There are two 
reasons for this. First, the eigenvalues of a nonsymmetric matrix can be very sensitive to 
small changes in the matrix elements. Second, the matrix itself can be defective, so that there 
is no complete set of eigenvectors. 

There are several basic building blocks in the QR algorithm, which is generally regarded 
as the most effective algorithm, for solving all eigenvalues of a real, unsymmetric matrix. 
These basic components of the QR algorithm are reviewed in Section II. Basic techniques 
to exploit the vector speed and multiple processors offered by modern high-performance 
computers are explained in Section III. An analysis of the Hessenberg reduction component 
in the QR algorithm is given in Section IV where both vector and parallel techniques are 
incorporated into the Hessenberg reduction component. Numerical examples are provided 
in Section V to evaluate the performance of the proposed method over the existing one. 
Conclusions and recommendations are given in Section VI. Finally, a listing of the 
Hessenberg reduction algorithm (in the form of Fortran coding) is provided in the appendix. 

II. Basic Components of the QR Algorithm [3,5] 

2.1 Balancing: 


1 


The idea of balancing is to use similarity transformations to make corresponding rows and 
columns of the matrix have comparable norms, thus reducing the overall norm of the matrix 
while leaving the eigenvalues unchanged. 

The time taken by the balanced procedure is insignificant as compared to the total time 
required to find the eigenvalues. For this reason, it is strongly recommended that a 
nonsymmetric matrix need to be balanced before even attempting to solve for eigen- 
solutions. 

2.2 Reduction to Hessenberg form: 

The strategy for finding the eigensolution of an unsymmetric matrix is similar to that of 
the symmetric case. First we reduce the matrix to a simpler Hessenberg form, and then we 
perform an iterative procedure on the Hessenberg matrix. An upper Hessenberg matrix has 
zeros everywhere below the diagonal except for the first subdiagonal. For example, in the 
6x6 case, the nonzero elements are: 

X X X X X X 

X X X X X X 

0 X X X X X 

0 0 X X X X 

0 0 0 X X X 

0 0 0 0 X X 


Thus, a procedure analogous to Gaussian elimination can be used to convert a general 
unsymmetric matrix to an upper Hessenberg matrix. The detailed coding of the Hessenberg 
reduction procedure is listed in subroutine OELMHS of the appendix. 

Once the unsymmetric matrix has already been converted into the Hessenberg form, the 
QR algorithm [3,5] itself can be applied on the Hessenberg matrix to find all the real and 
complex eigenvalues. For completeness, detailed coding of the QR algorithm on the 
Hessenberg matrix is listed in subroutine HQR of the appendix. 



In this section, a simple example of matrix times vector is used to explain some basic 
vector and parallel techniques which are useful for Hessenberg reduction algorithm. 

2-1 0 

Given a 3x3 Matrix A = -1 2-1 and a vector x = {1,0,0} r 

0 -1 1 


Here, the dimension of the system is N=3. The objectives are to develop efficient parallel - 
vector matrix times vector subroutines. 
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3.1 Row-by-Row conventional approach: 

DO 1 7=1,7/ 

DO 2 1=1^ 

B(I) = B(I) +A(I,J) *x(J) 

2 Continue 
1 Continue 


It should be emphasized here that in this approach, the value of B(I) corresponds to the 
final answer. 


3.2 Column-by-Column conventional approach: 

DO 1 J = 1 ,N 
DO 2 I = 1 ,N 
5(7) = B(I) + A(I,J) * x(J) 

2 Continue 
1 Continue 


It should be emphasized here that in this approach, the value of B(I) does NOT 
correspond to the final answer. B(I) only gives the partial (or incomplete) answer and it will 
give the final answer only if all values of J have been executed. It is also observed that x(J) 
is a constant (with respect to loop 2), thus the operations involved in loop 2 can be stated 
generally as: A new vector B = Old vector B + Constant * another vector A. 

3.3 Row-by-Row “vector unrolling” approach: 

Assuming the dimension N of the system is large, say N = 600, then the algorithm in 
Section 3.1 can be modified to improve the vector speed as following: 

NUNROL = 2 

DO 1 / = 1 ,N, NUNROL 

DO 2 J = 

5(7) = B(I) + A(I,J) * x(J) 

B(I+ 1) = 5(7+1 ) + A(/+1,J) * x(J) 

2 Continue 
1 Continue 


The operations involved inside loop 2 is referred to as “dot product” operations. 

3.4 Column-by-Column “loop-unrolling” approach 

The algorithm in Section 3.2 can be modified to improve the vector speed performance 
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NUNROL = 2 

DO 1 J = ^,N, NUNROL 

DO 2 I = 1,N 

B(I) = B(I) + A{I,J) * x(J ) 

+ /4(/J+1) * x(J + 1) 

2 Continue 
1 Continue 


The operations involved inside loop 2 is referred to as "saxpy" operations. 

3.5 Parallel-vector loop-unrolling approach: 

For multiple processors, the algorithm in Section 3.4 can be modified to take advantage 
of parallel speed (in addition to vector speed) 

NUNROL = 2 

Parallel DO 1 J = 1 ,N, NUNROL 

DO 2 I = 1 ,N 

B(I) = B(I) + A(I,J) * x (J) 

+ A(I,J+ 1) * x(J + 1) 

2 Continue 
1 Continue 

In this algorithm, each value of the index J (of loop 1) is assigned to different processors for 
parallel computation. 

IV. An Analysis of the Hessenberg Reduction Algorithm 

A careful look into the Hessenberg reduction algorithm of Section 2.2 and subroutine 
OELMHS of the appendix will reveal that the most intensive computations of Subroutine 
OELMHS occur in loops 140 and 150 of the code. Furthermore, the Fortran statement 
inside loop 150 can be generally expressed as: 

A(J, M) = A(J,M) + Y * A(J,I) 

or 

A new vector A(J, -) = old vector A (J, -) + (a constant) * another vector A(J,*) 

Thus, one can immediately see the similarity between loops 160 & 150 of Subroutine 
OELMHS and loops 1 & 2 of the matrix times vector algorithm presented in Section 3.2. 
From the experience we have had in section 3.5, we can therefore similarly apply the parallel 
computations in loop 160 and loop-unrolling (here NUNROL = 8 is used) for vector 
computations in loop 150 of subroutine OELMHS. 

For completeness, the entire parallel-vector version of the Hessenberg reduction, and the 
original QR algorithms are listed in the Appendix. 
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V. Numerical Examples 


In order to evaluate the numerical accuracy and the performance of the new parallel- 
vector Hessenberg Reduction portion of the QR algorithm, the following numerical tests are 
performed. 

Example 1: 

Find all eigenvalues of the following 2x2 unsymmetric matrix 

A = [2 -el 

[8 ij 

The analytical eigen-value solution for this problem is: 

A.., =1.5 + 6.91 i 
k 2 = 1.5 - 6.91 i 


which also matches with the computer solution. 

Example 2: 

In this example, the unsymmetric matrix [A] NxN is automatically generated for any dimension 
N of the matrix [A] (please refer to the code given in the Appendix). The accuracy and the 
performance of the new parallel-vector Hessenberg reduction algorithm is compared to the 
original subroutine. Since the QR algorithm itself is highly sequential, no attempts to 
parallelize and vectorize the QR algorithm have been made. However, the total solution 
time of the complete unsymmetric eigensolution process (= Hessenberg Reduction Time and 
QR Time) are also presented in Tables 1 and 2. 
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Table 1 : Vector Performance on the Alliant Using etime (t), fortran -DAS -O -alt -1 -OM 
where: 

- 1 option will tell which loop does not vectorize 

- OM option will not print warning messages 


“Original” CSI version 

HR - Hessen berg 
Reduction Time 
k QR Time , 


Size N 


“New” version 













Table 2 : Parallel- Vector Performance on Cray-YMP (Reynolds) Using tsecnd (). 


Size N 

“Original” CSI v 
(HR - Hessenberg' 
Reduction Time 
( QR Time 

1 Cray-YMP Process 

ersion 

or 

“New” version 

1 Cray-YMP 
Processor 

2 Cray-YMP 
Processors 

3 Cray-YMP 
Processors 

100 x 100 

[ 0.02 sec^ 
( 0.07 sec, 


[ 0.02 seel 
( 0.07 secj 

[ 0.03 sec'j 
( 0.07 secj 

[ 0.03 seel 
( 0.07 secj 

200 x 200 


r 0 . 12 ^ 

, 0 . 42 , 



'o.i -T 
, 0 . 42 , 



' 0 . 08 1 

, 0 . 41 ; 



' 0 . 0 / 

, 0 . 41 , 






1 


1 

1 


1 



1 

600 x 600 






1 



1 

1 


1 

800 x 800 




1 



1 


I 

1 




VI. Conclusions and Recommendations : 

The most popular and effective procedure to solve all eigenvalues of an unsymmetric 
matrix involved 2 major tasks, namely Hessenberg reduction form and QR algorithm on the 
Hessenberg matrix. In general, QR algorthm requires between 2 to 3 times more 
computational effort than the Hessenberg reduction algorithm. 

In this study, the parallel and vector speeds of the Hessenberg reduction algorithm has 
been developed and implemented on the Alliant and Cray-YMP (Reynolds) computers. 
Numerical results have indicated that the proposed parallel-vector Hessenberg reduction 
algorithm does offer computational advantages (without losing its accuracy) as compared to 
the existing algorithm. The time saving is more significant as the problem size increased. 
Further research work is critically needed to improve the unsymmetric eigensolution 
procedure (using the QR, or another better, new parallel algorithm). 
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APPENDIX 


Parallel-Vector Hessenberg Reduction And Sequential 
QR Algorithm [5] 
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FILE: UNSEIG FRC 


A1 OLD DOMINION UNiVERhl I Y 


C PARALLEL/VECTOR UNSYMMETRIC E I GENSOLVER by Qin & Nguyen, May 1992 ** 

c This is a working version of "unsymmetr i cal" eigen-solver 

c on the sun386 work station. On the Cray-YMP (Reynold or Sabre), 

c this "exact 11 same version should offer good vector & parallel 

c speed (only for subroutine to perform Hessenberg reduction). 

c..'..*For SMALL problems, the improvements due to par a 1 1 e 1 -vector 

c Hessenberg is NOT MUCH. However, for LARGE problems, since the 

c .Hessenberg reduction timing becomes more important (as compared to 

c the TOTAL eigen-solution time), the total time saving for the entire 

c ..... .e i gen-sol ut i on process is also very significant. 

c Since this version was developed specifically for CSI applications 

c (according to Peiman’s spec i f i cat i ons/requ i rements) , ALL EIGENVALUES 

c (and NONE of the corresponding EIGENVECTORS) of an N by N squared 

c unsymmatr i ca 1 matrix are found. 

c M ART I F I C I AL 11 datas of varous sizes (N ■ 2 > 800) with ALL REAL 

c and MIXED REAL & COMPLEX eigenvalues have been verified (by comparing 

c the results obtained by the original unsym. eigen-sol. taken from 

c ORACLE and the modified version from the ODU team, and also by HAND 

c CALCULATION for the size N = 2) 


Force PVQR of NP ident ME 
Shared REAL A (1000000) ,WK (1000,2) 

Shared REAL ER (1000) , E I (1000) , E I G (1000) 

Shared REAL EPS , ERRCK 

Shared INTEGER N , NM, NMM, NMAX , NST, MQ, I MODE , I ERR, nguyen 
End Declarations 

C *** THIS IS THE PROGRAM CALL UNSYMMETRIC E I GENSLVER ******** 
Barr i er 

WRITE (*,*) 'N, IM0DE (0=old vers i on) , nguyen (l=duc-s data) »■ 

READ (5**) N , i mode, nguyen 

WR I TE (*, *) 1 N IMODE NGUYEN = 1 , N , i MODE , nguyen 

ERRCK= 0.0000001 

eps=geteps ( i beta , i t , i rnd) 

wr i te (*,*) 1 *** EPS =‘,eps 

wr i te (*, 101) N, imode 

101 FORMAT (//,* INPUT PARAMETERS: 1 ,/, 

1 ’N = 1 , 1 5 » 1 ~ Size of System 1 //* 

1 1 I MODE = ’,15/ " = 0 is old sequential 1 //) 

End Barrier 

Forcecal 1 RESV (N , N , A , ER, E I , WK , I ERR, E I G , I MODE , nguyen) 

Joi n 
END 

FUNCTION GETEPS (IBETA, IT, IRND) 
a • 1.0 

10 a = a + a 

i f ( ( (a+1 .0) -a) -1 .0.eq.O.OO) go to 10 
b=l .0 

20 b=b+b 

if ( (a+b) -a.eq.0.00) go to 20 
q i na= (a+b) -a 
i beta=i nt (q i na) 
beta=f loat (ibeta) 
i t=0 
b= 1 .0 

30 i t= i t+1 

b=b*beta 
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i f ( ( (b+ 1 . 0) -b) - 1 . 0 . eq . 0 . 00) go to 30 
i rnd=0 

betaml=beta- 1 .0 
i f ( (a+betaml) -a.ne.O.OO) i rnd=l 
betai n=l .O/beta 
a=l .0 

do 40 i=l , i t+3 
a=a*beta i n 
40 continue 

50 if ((1 .0+a) -1 .O.ne.O.OO) go to 60 

a=a*beta 
go to 50 
60 eps=a 

i f ( (i beta . eq * 2) .or . (i rnd . eq .0) ) go to 70 

a=(a*(l . 0+a) ) / (1 .0+1 .0) 
if ((1 .0+a) -1 .0.ne.0.00) eps=a 
70 geteps=eps 

return 
end 

Forcesub RESV (MAX,N , A, ER, E I ,WK, I ERR, E I G, I MODE .nguyen) of NP 
$ ident ME 

INTEGER MAX.N, I ERR, IMODE 
Shared Integer LOW,IGH,NACC 
C F2.4 
c rtrtrtrt 

C FUNCTION - COMPUTES ALL THE EIGENVALUES AND SELECTED 

C PARAMETERS MAX - MAXIMUM ROW DIMENSION OF A 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


N - ORDER OF A 

A (MAX , N) - INPUT MATRIX (DESTROYED) 

ER (N) - CONTAINS REAL PART OF THE EIGENVALUES 

El (N) - CONTAINS IMAGINARY PART OF THE EIGENVALUES 

WK(-) - WORKING STORAGE OF FOLLOWING DIMENSION 

DIMENSION 3*N IF ISV+ILV = 0 

DIMENSION N*(N+7) OTHERWISE 
I ERR - INTEGER ERROR CODE 

= 0 NORMAL RETURN 

= -J J-TH EIGENVECTOR DID NOT CONVERGE. 

VECTOR SET TO ZERO. IF FAILURE OCCURS 
MORE THAN ONCE, INDEX FOR LAST 
OCCURRENCE IN I ERR. 

= J J-TH EIGENVALUE HAS NOT BEEN 

DETERMINED AFTER 30 ITERATIONS 

OUTPUT FORMAT - EIGENVALUES ARE STORED IN ASCENDING MAGNITUDE 

WITH COMPLEX CONJUGATES STORED WITH POSITIV 
IMAGINARY PARTS FIRST. THE EIGENVECTORS ARE 
PACKED AND STORED IN V IN THE SAME ORDER AS 
THEIR EIGENVALUES APPEAR IN ER AND El. 

ONLY ONE EIGENVECTOR IS COMPUTED FOR COMPLE 
CONJUGATES (FOR CONJUGATE WITH POSITIVE 
IMAGINARY PART). UPON ERROR EXIT -J, EIGEN- 
VALUES ARE CORRECT AND EIGENVECTORS 
ARE CORRECT FOR ALL NON-ZERO VECTORS. 

UPON ERROR EXIT J, EIGENVALUES ARE CORRECT 
BUT UNORDERED FOR INDICES I ERR+1 , I ERR+2 , . . . 
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C N AND NO EIGENVECTORS ARE COMPUTED. 

C REQUIRED ROUTINES - QXZl46,QXZl47,QXZ152 

C **** 

REAL A (N , N) , ER (N) , E I (N) , WK (N , 2) , E I G (1 ) 

End Declarations 

C DIMENSION A(MAX,N) ,ER(N) , E I (N) ,V(MAX,*) ,WK(N,*) 

c LOGICAL LTESTV 

c EQUIVALENCE (TESTV, LTESTV) 

CQIN DATA TRUE, FALSE / ' 77777777777777777777 1 o , 'OOOOOOOOOOOOOOOOOOOO ' 0 
CQIN +/ 

c DATA TRUE, FALSE / 7777777777777777777-0,0.0000000000000000000 / 

C 

C PRELIMINARY REDUCTION 

c **** 

Barr i er 
DO 2 J=1 , N 
DO 1 1*1, N 
i f (i . 1 1 * j) then 

a (i,j)=l. 373737373737/ (float (i+j)) 
el se 

A (i,j) *0,9731 9731973V (float (i+j+j/2)) 

end i f 


1 continue 

2 continue 
do 3 i“l » n 

3 a (i , I) afloat (i *i) 

c Due T. Nguyen added this portion to test '‘complex" eigen-solution ! 


i f (nguyen . eq . 1 ) then 
DO 29 J*1,N 
. DO 19 1-1, N 
i f (i . 1 1 . j) then 

a (i , j) — 1 .373737373737*10.0/ (float (i+j) ) 
else 

A (i ,j) =0.9731 9731 9731 * 10 . 0 / (float (i+j+j/2)) 
end i f 

19 continue 

29 continue 

do 39 i-1 »n 
39 a (i , i) =f loat (i) 


a (1 , 1) =2 . 
a (1 , 2) =-6 . 
a (2 , 1) =8 . 
a (2 , 2) =1 . 
end i f 


C ***** SAVE A FOR NORM CHECK ***** 

1 ow= 1 
i gh=n 
TIME0=0.0 
End Barrier 
tOO=TSECND 0 

c CALL QXZ146 (MAX, N, A, LOW, IGH.WK) 
tl 1=TSECND 0 
i f (imode.ne.O) then 
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1122 


Forcecal 1 QXZlW (MAX , N , LOW, I GH , A , WK (1 , 2) , e i g) 
el se 

Forcecal 1 OELMHS (MAX, N, LOW, I GH , A , WK (1 , 2) ) 
end i f 

T22=TSECND () 

T I ME0=T I MEO+T22-TOO 
write (6,*) '** ME CPU in 
wr i te (6,*) 1 ** ME CPU in 
i f (me.eq. 1) then 

wri te (*,*) '*** A 

do 1122 i=n-10,n 
wr i te (*,*) 1 A (' , i , ' ,n) = 
conti nue 
end i f 


QXZ 1 i+6 = 1 ,ME , 

QXZ147 (OELMHS) 


T1 1 -TOO 

= 1 , ME ,T22-T1 1 


>Wc>'c 


,a (i ,n) 


C 

C COMPUTE ALL EIGENVALUES AND NO EIGENVECTORS 
C **** 

Barr i er 
tOO=TSECND 0 
i f (i mode . eq .0) then 

call HQR (MAX, N, LOW, I GH , A , ER , E I , I ERR) 
else 

call qxz 152 1 (max, n , low, i gh , A , er , ei , i err) 
end i f 

tl 1=TSECND 0 

wr i te (*,*) 1 ** IMODE ,CPU time in QXZ152 = ' , imode, tl 1 -tOO 
i f (me.eq. 1) then 

write(*,*)' *** E i gen val ue#, real ER(I), imaginary E I ( I ) ***' 
do 7 I =n- 1 0 , n 
write(*,*) I ,er (i) ,ei (i) 

7 continue 

c rearrange eigenvalues according to ascending order (of real part) 

call ascend (n , er , e i ,wk) . 
end i f 

End Barrier 

RETURN 

END 

C — SUBPROGRAM QXZ 1 i+6 — FORMERLY KNOWN AS ROUTINE BALANC — 

C 

c 

SUBROUTINE QXZ 146 (NM, N , A , LOW, I GH , SCALE) 

C 

INTEGER I , J , K, L , M, N , J J , NM, I GH , LOW, I EXC 
REAL A (N , N) .SCALE (N) 

REAL C,F,G,R,S,B2, RADIX 
LOGICAL NOCONV 
C 

C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE BALANCE, 

C NUM. MATH. 13, 293-304(1969) BY PARLETT AND REINSCH. 

C HANDBOOK FOR AUTO. COMP., VOL . I I -L I NEAR ALGEBRA, 3 15“326 (1 97 1 ) - 
C 

C THIS SUBROUTINE BALANCES A REAL MATRIX AND ISOLATES 
C EIGENVALUES WHENEVER POSSIBLE. 

C 

C ON INPUT 
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NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL 
ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM 
DIMENSION STATEMENT. 

N IS THE ORDER OF THE MATRIX. 

A CONTAINS THE INPUT MATRIX TO BE BALANCED. 

ON OUTPUT 

A CONTAINS THE BALANCED MATRIX. 

LOW AND IGH ARE TWO INTEGERS SUCH THAT A ( I , J) 

IS EQUAL TO ZERO IF 

(1) I IS GREATER THAN J AND 

(2) J=1 LOW-1 OR I = I GH+1 N. 

SCALE CONTAINS INFORMATION DETERMINING THE 
PERMUTATIONS AND SCALING FACTORS USED. 

SUPPOSE THAT THE PRINCIPAL SUBMATRIX IN ROWS LOW THROUGH IGH 
HAS BEEN BALANCED, THAT P (J) DENOTES THE INDEX INTERCHANGED 
WITH J DURING THE PERMUTATION STEP, AND THAT THE ELEMENTS 
OF THE DIAGONAL MATRIX USED ARE DENOTED BY D(I,J). THEN 
SCALE (J) - P(J), FOR J = 1 ,..., LOW-1 
= D (J, J) , J = LOW, .... IGH 

* P (J) J = I GH+1,..., N. 

THE ORDER IN WHICH THE INTERCHANGES ARE MADE I S N TO I GH+1, 

THEN 1 TO LOW-1. 

NOTE THAT 1 IS RETURNED FOR IGH IF IGH IS ZERO FORMALLY. 

THE ALGOL PROCEDURE EXC CONTAINED IN BALANCE APPEARS IN 
QXZH6 IN LINE. (NOTE THAT THE ALGOL ROLES OF IDENTIFIERS 
K, L HAVE BEEN REVERSED.) 

C 

C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, 

C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY 

C 

C THIS VERSION DATED AUGUST 1 983 • 

C 

C BASED ON THE EISPACK VERSION 3 ROUTINE BALANC, AS MODIFIED 
C BY COMPUTER SCIENCES CORPORATION, MAY 1981*. 

C 

C — 

C 

RADIX = 16.00 
C 

B2 = RADIX * RADIX 
K = 1 
L = N 
GO TO 100 


C IN-LINE PROCEDURE FOR ROW AND 

C COLUMN EXCHANGE 
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20 SCALE (M) = J 

IF (J .EQ. M) GO TO 50 
C 

DO 30 I = 1, L 
F = A ( I , J) 

A (I ,J) = A (I ,M) 

A ( I , M) = F 
30 CONTINUE 
C 

DO 40 I = K, N 
F = A (J , I ) 

A (J , I ) = A(M, I) 

A (M, I ) = F 
1*0 CONTINUE 
C 

50 GO TO (80,130) , I EXC 


C SEARCH FOR ROWS ISOLATING AN EIGENVALUE 

C AND PUSH THEM DOWN 

80 I F (L .EQ. 1) GO TO 280 
L = L - 1 

C FOR J = L STEP -1 UNTIL 1 DO — 


100 DO 120 JJ = 1, L 
J = L + 1 - JJ 

DO 110 I = 1, L 

IF (I .EQ. J) GO TO 110 
IF (A (J , I ) .NE. 0.00) GO TO 120 
110 CONTINUE 

M = L 
I EXC = 1 
GO TO 20 
120 CONTINUE 

GO TO 140 


C SEARCH FOR COLUMNS ISOLATING AN EIGENVALUE 

C AND PUSH THEM LEFT 


130 K = K + 1 
C 

140 DO 170 J = K, L 
C 

DO 150 I = K, L 

IF (I .EQ. J) GO TO 150 
IF (A ( I , J) .NE. 0.00) GO TO 1 70 
150 CONTINUE 
C 

M = K 
I EXC = 2 
GO TO 20 
170 CONTINUE 


C NOW BALANCE THE SUBMATRIX IN ROWS K TO L 

DO 180 I = K, L 
180 SCALE (I) = 1.00 

C ITERATIVE LOOP FOR NORM REDUCTION 

190 NOCONV = .FALSE. 
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C 

DO 270 I = K, L 
C = 0.00 
R = 0.00 
C 

DO 200 J = K, L 

IF (J .EQ. I) GO TO 200 
C = C + ABS (A (J, I) ) 

R = R + ABS (A (I ,J)) 


200 CONTINUE 

C GUARD AGAINST ZERO C OR R DUE TO UNDERFLOW 

IF (C .EQ. 0.00 .OR. R .EQ. 0.00) GO TO 270 
G = R / RADIX 
F = 1 .00 
S = C + R 

210 IF (C .GE. G) GO TO 220 

F = F * RADIX 
C = C * B2 
GO TO 210 

220 G = R * RADIX 

230 IF (C .LT. G) GO TO 240 

F = F / RADIX 
C = C / B2 
GO TO 230 

C NOW BALANCE 

240 IF ( (C + R) / F .GE. 0.950 * S) GO TO 270 
G = 1 .00 / F 


SCALE (I) « SCALE (I) * F 
NOCONV = .TRUE. 

C 

DO 250 J = K, N 
250 A (I , J) = A (I , J) * G 
C 

DO 260 J = 1, L 
260 A (J , I ) = A(J, I) * F 
C 

270 CONTINUE 
C 

IF (NOCONV) GO TO 1 90 
C 

280 LOW = K 
IGH - L 
RETURN 

C ********** LAST CARD OF QXZ146 ********** 

END 

C — SUBPROGRAM QXZ147 — FORMERLY KNOWN AS ROUTINE ELMHES — 
C 

C 

Forcesub QXZ 1 47 (NM, N , LOW, I GH , A, I NT, temy) of NP ident ME 
C 

INTEGER N , NM, I GH , LOW, I NT (1) 

REAL A (N , N) , temy (1) 

Shared INTEGER LA , KP1 ,MM1 , MP1 , I AM 
Shared Logical ilock 
Shared REAL X,Y,XMUL,XMUL1 
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Private Real tema(lOOO) 

End Declarations 

C — 

c 

I AM=1 
Barrier 
LA - IGH - 1 
KP1 = LOW + 1 

IF (LA .LT. KP1) GO TO 200 
C 

End Barrier 
DO 180 -M = KP1 , LA 
Barrier 
End Barrier 
I F (ME . EQ. I AM) THEN 
MM1 = M - 1 
X = 0.00 
I = M 
C 

DO 100 J = M, IGH 

IF (ABS (A (J , MM1 ) ) .LE. ABS (X) ) GOTO 100 
X = A (J , MM 1 ) 

I = J 


100 CONTINUE 

C 

I NT (M) = I 

IF (I .EQ. M) GO TO 130 

C INTERCHANGE ROWS AND COLUMNS OF A 

DO 110 J = MM1 , N 
Y = A (I ,J) 

A (I , J) = A (M, J) 

A (M, J) = Y 
110 CONTINUE 


c 

DO 120 J = 1, IGH 
Y = A (J , I ) 

A (J , I ) = A (J , M) 
A (J , M) = Y 


120 CONTINUE 

C END INTERCHANGE . 

C130 IF (X .EQ. 0.00) GO TO 180 
130 CONTINUE 


END I F 
Barrier 
End Barrier 
Bar r i er 
i am=iam+l 

If (iam.gt.NP) i am=l 
End Barrier 

I F (X.EQ.O.OO) GO TO 1800 
I F (ME - EQ. I AM) THEN 
do 1301 i =m+l , i gh 
temy (i) =a (i ,mml) /x 

i f (temy (i) .ne.O.OO) a (i ,mml) =temy (i) 
c if (a (i ,mml) .eq.O.OO) then 
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c temy(i)=0.0 

c else 

c temy (i) =a (i ,mml) /x 

c a (i ,mml) = temy (i) 

c end if 

1301 continue 

C 

c DO 160 I = MP1 , IGH 

c Y = A (I ,MM1) 

c IF (Y .EQ. 0.00) GO TO 160 

c Y = Y / X 

c A (I , MM1 ) = Y 

C 

END I F 

do 1399 j-l.iflh 

1399 tema(j)=0.0 
jend= ( (i gh-m) /8) *8 

Barr i er 
I am= i am+1 

if (I AM.GT.NP) I AM=1 
End Barrier 
Barr i er 
End Barrier 

c do 1400 j j=m+l ,m+j end , 8 

Presched DO 1400 j j=m+l ,m+jend,8 

CD I R$ IVDEP 

do 1401 j=l,m-l 

c a (j , m) =a (j ,m) +temy (j j) *a (j , j j)+temy (jj+1) *a (j , j j + 1) 

tema (j) =tema (j)+temy (j j) *a (j , jj)+temy (j j + 1) *a (j , jj + 1) 

1 +temy (j j+2) *a (j , j j+2)+temy (j j+3) *a (j » jj+3) 

2 +temy (jj+4) *a (j , j j+4) +temy (jj+5) *a (j ,jj+5) 

3 +temy (jj+6) *a (j , j j+6)+temy (jj+7) *a (j , j j+7) 

1401 continue 

1400 End Presched DO 

Barr i er 

End Barrier 

Presched DO 1402 j j=jend+l+m, i gh 

CD I R$ IVDEP 

do 1403 j=l ,m-l 

1 403 tema (j)=tema (j)+temy (jj) *a (j , jj) 

1402 End Presched DO 

Barr i er 

End Barrier 
Critical i 1 ock 
do 14001 j=l ,m- 1 

14001 a(J,H) = a(J,rt) + tema(j) 

End Critical 
Barr i er 
End Barrier 

Presched DO 1411 jj=m+1,n 

CD I R$ IVDEP 

do 1412 i i =m+l , i gh 

H*12 a (i i , jj) =a (i i , jj) -a (m, jj) *temy (i i) 

1411 End Presched DO 

Barr i er 
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End Barrier 
IF (ME. EQ. I AM) THEN 
do 1 407 kk=m+1 , I gh 
xmul=temy (kk) 
xmul l=xmul *a (m, kk) 
c write (* , *) a (kk ,m) 

c a (kk ,m) =a (kk ,m) -temy (kk) *a (m,m) 

a (kk,m) =a (kk,m) -xmul*a (m,m) 

CD I R$ IVDEP 

do 1608 i k=m , kk 

1 608 a ( i k , m) =a ( i k , m) +xmu 1 *a ( i k , kk) 

CD I R$ IVDEP 

do 1 609 i k=kk+l , i gh 

1609 a (i k ,m) =a (i k,m) +xmul *a (i k , kk) +xmul l*temy (i k) 

c write (*,*) a (kk,m) , temy (kk) ,a (m,m) 

1407 continue 

END I F 

1800 continue 

c DO 140 J = M, N 

c 140 A (I ,J) = A (I , J) - Y * A (M, J) 

C 

c DO 150 J = 1, I GH 

c 150 A (J , M) = A (J ,M) + Y * A (J , I ) 

C 

c 160 CONTINUE 
C 

c801 write (*,*)' M-th step A(i,mml-n) igh.jend =' ,M, igh, jend 
c do 1500 I 1 = 1 , N 

c wr i te (*, 1501 ) (a (i i , j j) , j j=mml , n) 

c 500 continue 

c 501 format (lx, 10 (e9-3. lx) ) 

Barr i er 
I AM= I AM+1 

IF (I AM.GT.NP) I AM-1 
End Barrier 
180 CONTINUE 
C 

200 RETURN 

C ********** LAST CARD OF QXZ147 ********** 

C** THIS PROGRAM VALID ON FTN4 AND FTN5 ** 

END 

SUBROUTINE BALANC (NM,N,A,LOW, IGH, SCALE) 

C 

INTEGER I ,J,K,L,M,N,JJ,NM, IGH, LOW, IEXC 
REAL A (N , N) .SCALE (N) 

REAL C,F,G,R,S,B2, RADIX 
LOGICAL NOCONV 
C 

C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE BALANCE, 

C NUM. MATH. 1 3, 293"304 ( 1969 ) BY PARLETT AND REINSCH. 

C HANDBOOK FOR AUTO. COMP., VOL .1 I -L I NEAR ALGEBRA, 315-326 (1971). 

C 

C THIS SUBROUTINE BALANCES A REAL MATRIX AND ISOLATES 
C EIGENVALUES WHENEVER POSSIBLE. 

C 
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C ON INPUT 

C 

C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL 

C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM 

C DIMENSION STATEMENT. 

C 

C N IS THE ORDER OF THE MATRIX. 

C 

C A CONTAINS THE INPUT MATRIX TO BE BALANCED. 

C 

C ON OUTPUT 

C 

C A CONTAINS THE BALANCED MATRIX. 

LOW AND IGH ARE TWO INTEGERS SUCH THAT A (I , J) 

IS EQUAL TO ZERO IF 
C (1) I IS GREATER THAN J AND 

C (2) J=1 LOW-1 OR l-IGH+1, . .. ,N. 

C 

C SCALE CONTAINS INFORMATION DETERMINING THE 

C PERMUTATIONS AND SCALING FACTORS USED. 

C 

C SUPPOSE THAT THE PRINCIPAL SUBMATRIX IN ROWS LOW THROUGH IGH 
C HAS BEEN BALANCED, THAT P (J) DENOTES THE INDEX INTERCHANGED 
C WITH J DURING THE PERMUTATION STEP, AND THAT THE ELEMENTS 

C OF THE DIAGONAL MATRIX USED ARE DENOTED BY D(I,J). THEN 

C SCALE (J) = P(J), FOR J = 1,..., LOW-1 

C - D (J, J) , J = LOW, .... IGH 

C = P(J) J = IGH+1.....N. 

C THE ORDER IN WHICH THE INTERCHANGES ARE MADE IS N TO IGH+1, 

C THEN 1 TO LOW-1. 

C 

C NOTE THAT 1 IS RETURNED FOR IGH IF IGH IS ZERO FORMALLY. 

C 

C THE ALGOL PROCEDURE EXC CONTAINED IN BALANCE APPEARS IN 

C BALANC IN LINE. (NOTE THAT THE ALGOL ROLES OF IDENTIFIERS 

C K , L HAVE BEEN REVERSED.) 

C 

C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, 

C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY 

C 

C THIS VERSION DATED AUGUST 1983. 

C 

C 

C 

RADIX = 16.0E0 
C 

B2 = RADIX * RADIX 
K = 1 
L = N 
GO TO 100 


C IN-LINE PROCEDURE FOR ROW AND 

C COLUMN EXCHANGE 


20 SCALE (M) - J 

IF (J .EQ. M) GO TO 50 
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DO 30 I = 1, L 
F = A (I , J) 

A (I ,J) = A (I ,M) 

A ( I , M) = F 
30 CONTINUE 
C 

DO 40 I = K, N 
F = A (J , I ) 

A(J, I) = A (M, I) 

A (M, I ) = F 
40 CONTINUE 
C 

50 GO TO (80,130) , I EXC 


C SEARCH FOR ROWS ISOLATING AN EIGENVALUE 

C AND PUSH THEM DOWN 

80 IF (L .EQ. 1) GO TO 280 
L = L - 1 

C FOR J=L STEP -1 UNTIL 1 DO — 


100 DO 120 JJ = 1 , L 
J = L + 1 - JJ 

DO 110 I = 1, L 

IF (I .EQ. J) GO TO 110 
IF (A (J , I ) .NE. O.OEO) GO TO 120 
110 CONTINUE 

M = L 
I EXC = 1 
GO TO 20 
120 CONTINUE 

GO TO 140 


C SEARCH FOR COLUMNS ISOLATING AN EIGENVALUE 

C AND PUSH THEM LEFT 


130 K = K + 1 
C 

140 DO 170 J = K, L 
C 

DO 150 I = K, L 

IF (I .EQ. J) GO TO 150 
IF (A (I ,J) .NE. O.OEO) GO TO 170 
150 CONTINUE 
C 

M = K 
I EXC = 2 
GO TO 20 
170 CONTINUE 


C NOW BALANCE THE SUBMATRIX IN ROWS K TO L 

DO 180 I = K, L 
180 SCALE (I) = 1 .OEO 

C ITERATIVE LOOP FOR NORM REDUCTION 

190 NOCONV = .FALSE. 

C 

DO 270 I = K, L 
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C = O.OEO 
R = O.OEO 

DO 200 J = K, L 

IF (J .EQ. I) GO TO 200 
C = C + ABS (A (J, I)) 

R = R + ABS (A (I ,J)) 

200 CONTINUE 


C GUARD AGAINST ZERO C OR R DUE TO UNDERFLOW 

IF (C .EQ. O.OEO .OR. R .EQ. O.OEO) GO TO 270 
G = R / RADIX 
F = 1 .OEO 
S = C + R 

210 IF (C .GE. G) GO TO 220 
F = F * RADIX 
C = C * B2 
GO TO 210 

220 G = R * RADIX 
230 IF (C .LT. G) GO TO 21*0 
F = F / RADIX 
C = C / B2 
GO TO 230 

C NOW BALANCE 

2A0 IF ((C + R) / F .GE. 0.95E0 * S) GO TO 270 


G = 1 .OEO / F 

SCALE (I) - SCALE (I) * F 

NOCONV = .TRUE. 

C 


C 


C 

C 

C 


DO 250 J 
250 A (I , J) = 

DO 260 J 
260 A (J , I ) = 

270 CONTINUE 

IF (NOCONV) 


280 LOW = K 


= K, N 
A (I , J) * G 

- 1, L 
A (J, I) * F 


GO TO 190 


I GH = L 

RETURN 

END 

SUBROUTINE HQR (NM, N , LOW, IGH,H,WR,WI , I ERR) 

INTEGER I ,J,K,L,M,N,EN,LL,MM,NA,NM, I GH , ITN, ITS , L0W.MP2 , ENM2 , I ERR 
REAL H (N.N) ,WR (N) ,WI (N) 

REAL P,Q,R,S,T,W,X,Y,ZZ,N0RM,TST1 ,TST2 
LOGICAL NOTLAS 


THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE HQR, 
NUM. MATH. 1U, 219-231(1970) BY MARTIN, PETERS, AND WILKINSON. 
HANDBOOK FOR AUTO. COMP., VOL . I I -L I NEAR ALGEBRA, 359-3710971). 
C 

C THIS SUBROUTINE FINDS THE EIGENVALUES OF A REAL 
C UPPER HESSENBERG MATRIX BY THE QR METHOD, 
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c 

C ON INPUT 
C 

C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL 

C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM 

C DIMENSION STATEMENT. 

C 

C N IS THE ORDER OF THE MATRIX. 

C 

C LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING 

C SUBROUTINE BALANC. IF BALANC HAS NOT BEEN USED, 

C SET LOW=1 , I GH=N . 

C 

C H CONTAINS THE UPPER HESSENBERG MATRIX. INFORMATION ABOUT 

C THE TRANSFORMATIONS USED IN THE REDUCTION TO HESSENBERG 

C FORM BY ELMHES OR ORTHES, IF PERFORMED, IS STORED 

C IN THE REMAINING TRIANGLE UNDER THE HESSENBERG MATRIX. 

C 

C ON OUTPUT 
C 

C H HAS BEEN DESTROYED. THEREFORE, IT MUST BE SAVED 

C BEFORE CALLING HQR IF SUBSEQUENT CALCULATION AND 

C BACK TRANSFORMATION OF EIGENVECTORS IS TO BE PERFORMED. 

C 

C WR AND Wl CONTAIN THE REAL AND IMAGINARY PARTS, 

C RESPECTIVELY, OF THE EIGENVALUES. THE EIGENVALUES 

C ARE UNORDERED EXCEPT THAT COMPLEX CONJUGATE PAIRS 

C OF VALUES APPEAR CONSECUTIVELY WITH THE EIGENVALUE 

C HAVING THE POSITIVE IMAGINARY PART FIRST. IF AN 

C ERROR EXIT IS MADE, THE EIGENVALUES SHOULD BE CORRECT 

C FOR INDICES I ERR+1 , . . . , N . 

C 

C I ERR IS SET TO 

C ZERO FOR NORMAL RETURN, 

C J IF THE LIMIT OF 30*N ITERATIONS IS EXHAUSTED 

C WHILE THE J-TH EIGENVALUE IS BEING SOUGHT. 

C 

C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, 

C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY 

C 

C THIS VERSION DATED AUGUST 1983- 
C 

C 

c 

I ERR = 0 
NORM = O.OEO 
K = 1 

C STORE ROOTS ISOLATED BY BALANC 

C AND COMPUTE MATRIX NORM 

DO 50 I = 1, N 
C 

DO 40 J = K, N 

l»0 NORM *= NORM + ABS (H ( I , J) ) 

C 

K = I 
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IF (I .GE. LOW .AND. I .LE. I GH) GO TO 50 
WR (I) * H (I , I) 

Wl (I) = O.OEO 
50 CONTINUE 
C 

EN = I GH 
T = O.OEO 
ITN = 30*N 

C SEARCH FOR NEXT EIGENVALUES 

60 IF (EN .LT. LOW) GO TO 1001 
ITS = 0 
NA = EN - 1 
ENM2 = NA - 1 


C LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT 

C FOR L=EN STEP -1 UNTIL LOW DO — 


70 DO 80 LL - LOW, EN 

L = EN + LOW - LL 

IF (L .EQ. LOW) GO TO 100 

S = ABS (H (L-l , L- 1) ) + ABS (H (L , L) ) 

IF (S .EQ. O.OEO) S = NORM 
TST1 = S 

TST2 = TST1 + ABS (H (L , L- 1) ) 

IF (TST2 .EQ. TST1) GO TO 100 
80 CONTINUE 

C FORM SHIFT 

100 X = H (EN , EN) 

IF (L .EQ. EN) GO TO 270 
Y = H (NA , NA) 

W = H (EN , NA) * H (NA, EN) 

IF (L .EQ. NA) GO TO 280 
IF (ITN .EQ. 0) GO TO 1000 

IF (ITS .NE. 10 .AND. ITS .NE. 20) GO TO 130 

C FORM EXCEPTIONAL SHIFT 

T = T + X 
C 

DO 120 I = LOW, EN 
120 H ( I , I ) = H ( I , I ) -X 


C 

S = ABS (H (EN,NA)) + ABS (H (NA , ENM2) ) 

X = 0.75E0 * S 
Y = X 

W = -0.4375E0 * S * S 
130 ITS = ITS + 1 
ITN = ITN - 1 

C LOOK FOR TWO CONSECUTIVE SMALL 

C SUB-DIAGONAL ELEMENTS. 

C FOR M=EN-2 STEP -1 UNTIL L DO — 


DO 140 MM = L, ENM2 
M = ENM2 + L - MM 
ZZ = H (M , M) 

R = X - ZZ 
S = Y - ZZ 

P = (R * S - W) / H (M+ 1 , M) + H (M , M+ 1 ) 
Q = H (M+l , M+l) - ZZ - R - S 
R = H (M+2.M+1) 
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S = ABS(P) + ABS(Q) + ABS(R) 

P = P / S 
Q = Q / S 
R = R / S 

IF (M .EQ. L) GO TO 150 

TST 1 = ABS (P) * (ABS (H (M-l ,M-1) ) + ABS (ZZ) + ABS (H (M+l ,M+1) ) ) 
TST2 = TST1 + ABS (H (M,M-1) ) * (ABS (Q) + ABS (R) ) 

IF (TST2 .EQ. TST1) GO TO 150 
140 CONTINUE 
C 

150 MP2 = M + 2 
C 

DO 160 I = MP2, EN 
H (I , I -2) = O.OEO 
IF (I .EQ. MP2) GO TO 160 
H (I ,1-3) = O.OEO 
160 CONTINUE 


C DOUBLE QR STEP INVOLVING ROWS L TO EN AND 

C COLUMNS M TO EN 


DO 260 K = M, NA 

NOTLAS = K .NE. NA 
IF (K .EQ. M) GO TO 170 
P = H(K,K-1) 

Q = H (K+l , K- 1 ) 

R = O.OEO 

IF (NOTLAS) R = H(K+2,K-1) 

X = ABS (P) + ABS (Q) + ABS (R) 

IF (X .EQ. O.OEO) GO TO 260 
P = P / X 

Q - Q / X 

R = R / X 

170 S = SIGN (SQRT (P*P+Q*Q+R*R) ,P) 

IF (K .EQ. M) GO TO 180 
H (K, K- 1) = -S * X 
GO TO 190 

180 IF (L .NE. M) H (K, K- 1) - -H(K,K-1) 
190 P = P + S 
X = P / S 
Y = Q / S 
ZZ = R / S 

Q = Q / P 
R = R / P 


IF (NOTLAS) GO TO 225 

C ROW MODIFICATION 

DO 200 J = K, N 

P = H (K, J) + Q * H (K+l , J) 
H(K,J) = H(K,J) - P * X 
H (K+l , J) = H (K+l , J) - P * Y 
200 CONTINUE 
C 

J « Ml NO (EN , K+3) 

C COLUMN MODIFICATION .... 

DO 210 I - 1, J 

P - X * H(I,K) + Y * H (I ,K+1) 
H(I.K) = H ( I , K) - P 
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H ( I , K+ 1 ) = H ( I ,K+1) - P * Q 
210 CONTINUE 

GO TO 255 
225 CONTINUE 

C ROW MODIFICATION 


DO 230 J = K, N 

P = H(K,J) + Q * H (K+l , J) + R * H (K+2 , J) 
H (K, J) = H (K, J) - P * X 
H (K+l , J) = H (K+l , J) - P * Y 
H (K+2 , J) = H (K+2 , J) - P * ZZ 


230 CONTINUE 
C 

J = Ml NO (EN , K+3) 

C COLUMN MODIFICATION 


DO 240 I - 1, J 

P = X * H ( I , K) + Y * H (I , K+ 1 ) + ZZ * H ( I , K+2) 
H(I,K) = H ( I , K) - P 
H (I , K+l ) = H (I , K+l) - P A Q 
H (I , K+2) = H (I , K+2) - P * R 
240 CONTINUE 

255 CONTINUE 

C 

260 CONTINUE 
C 

GO TO 70 

C ONE ROOT FOUND 

270 WR (EN) = X + T 
Wl (EN) = O.OEO 
EN = NA 
GO TO 60 

C TWO ROOTS FOUND 

280 P = (Y - X) / 2.0E0 
Q = P * P + W 
ZZ = SQRT (ABS (Q) ) 

X = X + T 

IF (Q .LT. O.OEO) GO TO 320 

C REAL PAIR 

ZZ = P + S I GN (ZZ, P) 

WR (NA) = X + ZZ 
WR (EN) * WR (NA) 

IF (ZZ .NE. O.OEO) WR (EN) = X - W / ZZ 
Wl (NA) = O.OEO 
Wl (EN) = O.OEO 
GO TO 330 

C COMPLEX PAIR 


320 

WR 

(NA) - 

X + P 


WR 

(EN) - 

X + p 


Wl 

(NA) = 

ZZ 


Wl 

(EN) = 

-ZZ 

330 

EN 

= ENM2 


GO TO 60 


C SET ERROR — ALL EIGENVALUES HAVE NOT 

C CONVERGED AFTER 30*N ITERATIONS 


1000 I ERR = EN 

1001 RETURN 
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END 

Forcesub OELMHS (NM,N, LOW, I GH, A, RINDEX) of NP ident ME 
C 

REAL 

+ A (N,N) , RINDEX (I GH)' 

INTEGER 

+ I GH , LOW, N, NM 
REAL 
+ X, Y 

Shared INTEGER KP 1 , LA, MMI , MP 1 
End Declarations 
C 

C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE ELMHES, 

C NUM. MATH. 12, 3^9-368 ( 1 968 ) BY MARTIN AND WILKINSON. 

C HANDBOOK FOR AUTO. COMP., VOL . I I -L I NEAR ALGEBRA, 339-358 ( 1971 ). 

C 

C GIVEN A REAL GENERAL MATRIX, THIS SUBROUTINE 

C REDUCES A SUBMATRIX SITUATED IN ROWS AND COLUMNS 

C LOW THROUGH I GH TO UPPER HESSENBERG FORM BY 
C STABILIZED ELEMENTARY SIMILARITY TRANSFORMATIONS. 

C 

C ON INPUT 

C 

C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL 

C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM 

C DIMENSION STATEMENT. 

C 

C N IS THE ORDER OF THE MATRIX. 

C 

C LOW AND I GH ARE INTEGERS DETERMINED BY THE BALANCING 

C SUBROUTINE BALANC. IF BALANC HAS NOT BEEN USED, 

C SET LOW=l , I GH=N . 

C 

C A CONTAINS THE INPUT MATRIX. 

C 

C ON OUTPUT 

C 

C A CONTAINS THE HESSENBERG MATRIX. THE MULTIPLIERS 

C WHICH WERE USED IN THE REDUCTION ARE STORED IN THE 

C REMAINING TRIANGLE UNDER THE HESSENBERG MATRIX. 

C 

C RINDEX CONTAINS INFORMATION ON THE ROWS AND COLUMNS 

C INTERCHANGED IN THE REDUCTION. 

C ONLY ELEMENTS LOW THROUGH I GH ARE USED. 

C 

C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, 

C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY 

C 

C THIS VERSION DATED AUGUST 1983 - 

C 

c 

C 

LA - I GH - I 

KP 1 = LOW + 1 

IF (LA .LT. KP1) RETURN 
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C 


C 


100 

C 


C 


1 10 

c 


120 

C 

130 

C 


C 

140 

C 

150 

c 

160 

c 

180 

c 


c 
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Barr i er 

DO 180 M = KP 1 , LA 
MM1 = M - 1 
X = O.ODO 
I = M 

DO 100 J = M, IGH 

IF (ABS (A (J,MM1) ) .LE. ABS (X) ) GOTO 100 
X = A ( J , MM 1 ) 

I = J 
CONTINUE 

RINDEX(M) = REAL ( I ) 

IF (I .EQ. M) GO TO 130 

INTERCHANGE ROWS AND COLUMNS OF A 

DO 110 J = MM1, N 

Y = A (I , J) 

A ( I , J) = A (M, J) 

A (M, J) = Y 
CONTINUE 

DO 120 J = 1, IGH 

Y = A (J , I ) 

A (J , I ) = A(J,M) 

A (J , M) = Y 
CONTINUE 

END INTERCHANGE 

IF (X .EQ. O.ODO) GO TO 180 
MP1 = M + 1 

DO 1 60 I = MP1, IGH 

Y = A (I , MM1) 

IF (Y .EQ. O.ODO) GO TO 160 

Y = Y / X 

A ( I , MM1) = Y 

DO 140 J = M, N 
A (I , J) = A (I , J) - Y * A (M, J) 

DO 150 J = 1, IGH 
A (J ,M) = A (J ,M) + Y * A (J, I) 

CONTINUE 

CONTINUE 

End Barrier 
RETURN 
END 

SUBROUTINE QXZ1521 (NM, N , LOW, I GH . H,WR, Wl , I ERR) 

INTEGER I , J , K , L ,M, N , EN , LL , MM, NA , NM, IGH, ITN , ITS,L0W,MP2,ENM2, I ERR 
REAL H (N,N) ,WR(N) ,WI (N) 

REAL P,Q,R,S,T,W,X,Y,ZZ, N0RM.TST1 ,TST2 
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LOGICAL NOTLAS 
C 

I ERR = 0 
NORM = 0.00 
K = 1 

C STORE ROOTS ISOLATED BY QXZ146 

C AND COMPUTE MATRIX NORM 

DO 50 I = 1 , N 
C 

DO 40 J = K, N 

40 NORM = NORM + ABS (H (I , J) ) 

C 

K = I 

IF (I .GE. LOW .AND. I .LE. I GH) GO TO 50 
WR (I ) = H (I , I) 

Wl (I) = 0.00 
50 CONTINUE 
C 

EN = I GH 
T = 0.00 
ITN = 30*N 

C SEARCH FOR NEXT EIGENVALUES 

60 IF (EN .LT. LOW) GO TO 1001 
ITS = 0 
NA = EN - 1 
ENM2 = NA - 1 


C LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT 

C FOR L=EN STEP -1 UNTIL LOW DO — 


70 DO 80 LL = LOW, EN 

L = EN + LOW - LL 

IF (L .EQ. LOW) GO TO 100 

S = ABS (H (L-l ,L-1)) + ABS (H (L , L) ) 

IF (S .EQ. 0.00) S = NORM 
TST1 = S 

TST2 = TST1 + ABS (H (L , L- 1 ) ) 

IF (TST2 .EQ. TST1) GO TO 100 
80 CONTINUE 

C FORM SHIFT 

100 X = H (EN , EN) 

IF (L .EQ. EN) GO TO 270 

Y = H (NA, NA) 

W = H (EN , NA) * H (NA , EN) 

IF (L .EQ. NA) GO TO 280 
IF (ITN .EQ. 0) GO TO 1000 

IF (ITS .NE. 10 .AND. ITS .NE. 20) GO TO 130 

C FORM EXCEPTIONAL SHIFT 

wr i te (*, *) '** EN, T X =' ,EN,T,X 
T = T + X 
C 

DO 120 I = LOW, EN 
120 H (I , I) = H (I , I) - X 
C 

S = ABS (H (EN , NA) ) + ABS (H (NA, ENM2) ) 

X =0.750 * S 

Y = X 


29 


FILE: UNSEIG FRC 


ftl 0 L L) UOMINIUN UNIVtKbllY 


W = -0.43750 * S * S 
130 ITS = ITS + 1 
ITN = ITN - 1 

C LOOK FOR TWO CONSECUTIVE SMALL 

C SUB-DIAGONAL ELEMENTS. 

C FOR M=EN-2 STEP -1 UNTIL L DO — 


DO HO MM = L, ENM2 
M = ENM2 + L - MM 
ZZ = H(M,M) 

R = X - ZZ 
S = Y - 11 

P = (R * S - W) / H (M+ 1 , M) + H(M,M+1) 

Q = H(M+1,M+1) - 11 - R - S 
R = H (M+2.M+1) 

S = ABS (P) + ABS (Q) + ABS (R) 

P = P / S 
Q = Q / S 
R = R / S 

IF (M .EQ. L) GO TO 1 50 

TST1 = ABS (P) * (ABS (H (M— 1 ,M-1) ) + ABS (ZZ) + ABS (H (M+l ,M+1) ) ) 
TST2 = TST1 + ABS (H (M,M-1) ) * (ABS (Q) + ABS (R) ) 

IF (TST2 .EQ. TST1) GO TO 1 50 
140 CONTINUE 
C 

150 MP2 - M + 2 
C 

DO 160 I = MP2 , EN 
H (I , 1-2) = 0.00 
IF (I .EQ. MP2) GO TO 160 
H (I , I -3) =0.00 
160 CONTINUE 


C DOUBLE QR STEP INVOLVING ROWS L TO EN AND 

C COLUMNS M TO EN 


DO 260 K = M, NA 

NOTLAS = K .NE. NA 
IF (K .EQ. M) GO TO 170 
P = H (K , K- 1 ) 

Q = H (K+l , K — 1 ) 

R = 0.00 

I F (NOTLAS) R = H (K+2.K-1) 

X = ABS (P) + ABS (Q) + ABS (R) 

IF (X .EQ. 0.00) GO TO 260 
P = P / X 
Q = Q / X 
R = R / X 

170 S = S I GN (SQRT (P*P+Q*Q+R*R) ,P) 

IF (K .EQ. M) GO TO 180 
H(K.K-l) = -S * X 
GO TO 190 

180 IF (L .NE. M) H(K.K-l) = -H(K,K-1) 

190 P = P + S 
X = P / S 
Y = Q / S 
ZZ = R / S 
Q = Q / P 
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R = R / P 

IF (NOTLAS) GO TO 225 

C ROW MODIFICATION 

DO 200 J = K, N 

P = H(K,J) + Q' * H (K+l , J) 
H(K,J) = H (K , J) - P * X 
H (K+l , J) = H (K+l , J) - P * Y 
200 CONTINUE 
C 

J = Ml NO (EN,K+3) 

C COLUMN MODIFICATION .. 


DO 210 I = 1, J 

P - X * H (I ,K) + Y * H (I , K+l) 
H (I ,K) = H (I ,K) - P 
H (I , K+l) = H (I , K+l) - P * Q 


210 CONTINUE 

GO TO 255 
225 CONTINUE 

C ROW MODIFICATION 


DO 230 J = K, N 

P - H(K,J) + Q * H (K+l , J) + R * H (K+2, J) 
H (K , J) = H (K , J) - P * X 
H (K+l , J) = H (K+l , J) - P * Y 
H (K+2 , J) = H (K+2 , J) - P * ZZ 


230 CONTINUE 
C 

J = M I NO (EN , K+3) 

C COLUMN MODIFICATION 


DO 240 I = 1 , J 

P = X * H ( I , K) + Y * H (I .K+l) + ZZ * H ( I , K+2) 
H ( I , K) = H ( I , K) - P 
H (I , K+l) = H (I , K+l) - P * Q 
H (I , K+2) = H (I , K+2) - P * R 
2^0 CONTINUE 

255 CONTINUE 

C 

c wr i te (*, *) ' NOTLAS , K , H (K, K) = 1 , NOTLAS, K, H (K,K) 

260 CONTINUE 
C 

GO TO 70 

C ONE ROOT FOUND 

270 WR (EN) = X + T 
Wl (EN) = 0.00 
EN = NA 
GO TO 60 

C TWO ROOTS FOUND 

280 P = (Y - X) / 2.00 
Q = P * P + W 
ZZ = SQRT (ABS (Q) ) 

X = X + T 


c ***** the following if is added by Qin ** 

I F ( Q. LT.O.OO) go to 320 
C REAL PAIR 


ZZ - P + S I GN (ZZ, P) 
WR (NA) = X + ZZ 
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WR(EN) = WR(NA) 

IF (ZZ .NE. O.OO) WR(EN) = X - W / ZZ 
Wl (NA) = 0.00 
Wl (EN) = 0.00 
GO TO 330 

C COMPLEX PAIR 

320 WR (NA) = X + P 
WR (EN) = X + P 
Wl (NA) = ZZ 
Wl (EN) = -ZZ 
330 EN = ENM2 
GO TO 60 


C SET ERROR — ALL EIGENVALUES HAVE NOT 

C CONVERGED AFTER 30*N ITERATIONS 


1000 I ERR = EN 

1001 RETURN 
END 

subroutine ascend (n, er , e i ,wk) 
c implicit rea 1 *8 (a-h , o-z) 

real er (1) ,ei (1) ,wk (n, 1) 


do 1 i=l ,n 

smal ^999999999- 

do 2 j=l , n 

if( er (j) . 1 t.smal 1 ) then 
smal l=er (j) 

1 ocate=j 
wk (i , 1) =er (j) 
wk (i , 2) =e i (j) 
end i f 

2 continue 

er (locate) =999999999 • 

I continue 
c ..... . 

do 21 i=l ,n 
er (i) =wk (i , 1) 

21 ei (i)-wk(i ,2) 

c 

write(6,*) 'real 6 imaginary evalues in ascending order 1 

do 11 i=n-10,n 

write (6,*) i , er ( i ) , e i ( i ) 

II continue 
return 
end 
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